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Abstract 


A  finite  element  program  is  presented  which  computes 
displacements,  strains,  and  stresses  in  wood  members  of 
arbitrary  shape  which  are  subjected  to  plane  strain/stress- 
loading  conditions.  This  report  extends  a  program 
developed  by  R.  L.  Taylor  in  1977,  by  adding  both  the  cubic 
isoparametric  finite  element  and  the  capability  to  analyze 
nonisotropic  materials.  The  computer  subroutines  developed 
by  the  author  are  listed  in  this  report,  along  with  both  the 
details  for  incorporating  them  into  Taylor  s  program  and  the 
required  user  instructions. 

Keywords:  Finite  element  analysis,  computer  program, 
isoparametric  elements,  stress  analysis,  orthotropic 
materials,  anisotropic  materials,  plane  loading,  design,  cubic 
finite  element,  quadratic  finite  element. 
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Plane  Stress  Analysis 
of  Wood  Members 
Using  Isoparametric 
Finite  Elements 

A  Computer  Program 


By 

TERRY  D.  GERHARDT,  Research  Engineer 


Introduction 

The  finite  element  (FE)  computer  program  written  by  the 
author  and  presented  in  this  report  was  developed  as  part 
of  a  cooperative  research  effort  involving  the  National 
Wooden  Pallet  and  Container  Association  (NWPCA), 

Virginia  Polytechnic  Institute  ^nd  State  University  (VPI&SU), 
-aneHthe^ttSDA  F0resTT5ervtee^This  research  is  designed  to 
establish  rational  engineering  design  procedures  for  wood 
pallets.  The  author’ s  role  in  this  endeavor  is  to  determine 
the  stiffness  and  strength  of  notched  stringer  members  of 
pallets  as  functions  of  notch  geometry,  material  properties, 
and  loading  conditions.  As  part  of  this  effort,  the  author 
developed  the  FE  program  described  in  this  paper  to 
compute  displacements  and  stresses  in  wood  members  of 
any  geometrical  shape  which  are  under  arbitrary  plane 
stress  or  strain-loading  conditions.  This  computer  program 
is  to  be  applied  to  the  notched-stringer  problem.  Details  of 
the  program  development,  user  instructions,  and  program 
listing  are  presented  in  this  report.  The  program  was 
verified  by  a  comparison  of  FE  predictions  of  displacements 
and  strains  in  center-loaded,  double-tapered  wood  beams 
with  data  available  in  the  literature  This  comparison  is 
presented  in  another  paper  ^ — 

The  developed  subroutines  are  listed  in  the  appendix, 
however  no  other  support  is  offered. 


*  Maintained  at  Madtoon,  Wto.,  in  cooperation  with  the  University  of 
Wisconsin. 

*  Italicized  numbers  in  parentheses  refer  to  literature  cited  at  end  of  report. 

*  Gerhardt,  T.  D.  On  flnftt  aiement  modeing  of  tapered  wood  beams.  In 
preparation.  U.S.  Oep.  Agrtc.  Forest  Sec.,  For.  Prod.  Lab.,  Madison,  Wis. 


Program  Development 

The  decision  to  develop  an  FE  program  rather  than 
purchase  one  of  the  many  multipurpose  programs  available 
was  based  on  several  factors: 

(1)  The  existence  of  subroutines  in  the  literature  to  form 
the  basis  of  a  general-purpose  FE  program, 

(2)  The  belief  that  a  developed  program  could  be  more 
readily  expanded  for  future  research  needs,  and 

(3)  The  desire  to  include  the  cubic  isoparametric  element 
in  the  program. 

An  FE  computer  program  developed  by  Taylor  (6)  was  used 
as  a  starting  point.  Taylor  s  program  can  input  data  for 
one-,  two-,  or  three-dimensional  structures.  It  is  written  in  a 
modular  form:  Adding  a  new  element  requires  writing  a 
single  subroutine.  This  flexibility  is  quite  appealing  to  the 
researcher  because  the  added  element,  whether  a  three- 
dimensional  solid  element  or  even  a  fluid  or  heat-transfer 
element,  utilizes  existing  (and  debugged)  code  for  data 
input,  matrix  assembly,  matrix  inversion,  etc.  The  program 
has  a  macro  instruction  language  that  allows  a  variable 
algorithm  capability.  Also,  storage  requirements  are 
dynamically  assigned  for  each  problem  and  stored  in  a 
single  array.  In  this  manner,  available  computer  memory  is 
used  efficiently  for  any  type  of  problem.  Finally,  the 
following  isoparametric  plane  elements  are  available: 
triangular,  linear  quadrilateral,  quadratic  serendipity  and 
Lagrangian  quadrilateral  elements.  A  much  more  detailed 
description  of  the  Taylor  FE  program  is  given  in  (7). 


The  program  presented  here  adds  to  the  Taylor  program 
the  capability  to  analyze  nonisotropic  materials  and  the 
cubic  quadrilateral  isoparametric  plane  element.  The  former 
addition  makes  analysis  of  wood-  or  composite-based 
structures  possible  by  providing  proper  formulation  for 
elements  with  orthotropic  and  anisotropic  (wood  with  slope 
of  grain)  elastic  properties.  The  cubic  element  can  be  easily 
collapsed  to  provide  an  accurate,  fully  compatible  fracture 
mechanics  element  (5).  Although  the  cubic  element  is  not 
included  in  commercially  available  FE  programs,  some 
researchers  have  used  the  cubic  element  to  model  certain 
regions  in  notch  problems  (1,2). 

User  Instructions 

To  properly  input  data  for  an  FE  run,  first  review  Taylor’s 
(6)  instructions  in  section  24.3,  pages  690  to  695.  These 
instructions  need  to  be  modified  only  minimally,  with  the 
exception  of  inputting  the  material  property  data.  The 
following  comments  will  consider  both  material  property 
input  and  use  of  the  cubic  isoparametric  element. 

Inputting  Material  Property  Data. 

Four  cards  are  now  required  for  each  material  in  the 
structure  instead  of  the  three  required  by  Taylor  (6).  The 
first  two  cards  are  identical  to  Taylor’s.  The  first  card 
indicates  that  material  property  data  follows,  and  the 
second  card  inputs  the  material  set  number  and  the  element 
type  (I  EL).  I  EL  should  be  input  as  1  for  any  of  the  plane 
stress/strain  elements.  Card  3  is  in  a  4I5.F10.0  format  as 
follows: 

Columns  Variable 

1-5  MATYPE 

6-10  I 

11-15  L 

16-20  K 

21-30  ANGLE 

MATYPE  is  the  material-type  variable  and  has  the  following 
values: 

1  for  isotropic  materials. 

2  for  materials  orthotropic  in  the  global  x-y  plane. 

3  for  materials  orthotropic  in  a  local  1-2  plane  (see  fig.  1). 

4  for  anisotropic  materials. 

I  is  the  plane  loading  variable  and  has  the  following  values: 

*  0  for  plane  strain  loading. 

*  0  for  plane  stress  loading. 

L  is  the  order  of  Gaussian  quadrature  specified  for  stiffness 
matrix  determination  (L  x  L  points/element). 

K  is  the  order  of  Gaussian  quadrature  specified  for  stress 
determination  (K  x  K  points/element). 

ANGLE  is  the  counterclockwise  orientation  of  1-2  local 
coordinates  from  global  x-y  coordinates  (only  used  when 
MATYPE  *  3);  <?  in  figure  1. 

The  recommended  order  of  Gaussian  quadrature  (7)  is  L  * 

K  **  3  for  the  cubic  12-node  element,  L  *  K  *  2  for  the 
quadratic  8-node  element,  and  L  -  K  -  1  for  both  the 
linear  4-node  and  triangular  3-node  elements. 


Figure  1  — Principal  material  axes.  (Ml 51918) 

The  fourth  card  is  in  an  8F10.0  format.  The  form  of  the 
fourth  card  depends  on  the  value  of  MATYPE  specified  on 
the  third  card- 


MATYPE  =  1 

Columns 

Variable 

1-10 

D(4) 

11-20 

E 

21-30 

MATYPE  -  2 

V 

Columns 

Variable 

1-10 

D(4) 

11-20 

tr 

*“X 

21-30 

31-40 

E, 

41-50 

G«y 

51-60*4 

E,* 

61-70* 

»'«* 

71-80* 

MATYPE  =  3 

'V* 

Columns 

Variable 

1-10 

D(4) 

11-20 

E, 

21-30 

31-40 

e2 

41-50 

g,2 

51-60* 

61-70* 

^13* 

71-80* 

MATYPE  -  4 

*'23# 

Columns 

Variable 

1-10 

D(4) 

11-20 

D„ 

21-30 

d,2 

31-40 

D„ 

41-50 

^22 

51-60 

oa 

61-70 

^33 

D(4)  is  material  density— not  needed  for  static  analysis. 

E  is  modulus  of  elasticity. 

G  is  the  shear  modulus. 
v  is  Poisson’s  ratio. 

D,  are  componets  of  the  symmetric  ’Moduli  of  Elasticity” 
matrix  for  an  anisotropic  material.  These  components  are 
defined  in  (3). 


*  (*)  indicates  that  these  properties  ere  only  required  for  plane  strain  analysis 
(I  -  0). 


r 


Specif  tying  12-Node  Cubic  Isoparametric  Element 

The  following  steps  are  needed  to  specify  the  12-node  cubic 
isoparametric  element: 

(1)  In  the  control  card  set  NEN  —  12. 

(2)  When  inputting  element  data  (ELEM)  use  the  local 
node  numbering  sequence  shown  in  figure  2. 

The  cubic  element  can  be  degenerated  to  either  quadratic 
or  linear  displacement  fields  on  any  desired  side  by 
specifying  IX  ~  0  for  one  or  two  of  the  side  nodes.  This 
capability  is  useful  when  the  structure  is  to  be  modeled  by 
more  than  one  type  of  element.  For  example,  it  may  be 
desirable  to  model  regions  in  areas  of  suspected  stress 
concentrations  with  cubic  elements  and  the  remainder  of 
the  structure  with  quadratic  or  linear  elements. 

Implementation  of  Program  and  Listing 

The  details  required  for  incorporating  the  developed 
program  into  Taylor’s  program  (6)  are  described  in  this 
section. 

Deletions 

Delete  statements  31-36  in  subroutine  SHAP2  in  (6).  This 
eliminates  the  nine-node  Lagrangian  element. 

Additions 

Add  subroutines  ELPROP,  SHAP3,  TRANSF,  and  BTREB 
in  the  appendix  to  the  program  in  (6). 

Add  the  following  Common’  statement  to  subroutine 
PMESH  in  (6): 

COMMON  /  WRITE  /  PRT 

Add  the  following  ’Implicit’  statement  to  all  subroutines  in 
(6)  if  double  precision  computations  are  desired: 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

Substitutions 

Substitute  subroutine  ELMT01  in  the  appendix  for 
subroutine  ELMT01  in  (6). 


{-/,/}  a.n 

4  II  7  3 


HrD  drD 

({,  y)  SPACE 

Figure  2 — The  cubic  isoparametric  element  (local 
node  numbering).  (Ml 51 91 7) 
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Substitute  the  two  statements  that  follow  for  statement  18 
in  subroutine  SHAPE  in  (6): 

120  IF(NEL.GT.4.AND.NEL.LE.8)  CALL  SHAP2 
(SS,TT,SHP,IX,NEL) 

IF(NEL.GT.8)  CALL  SHAP3(SS,TT,SHP,IX) 
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Appendix 

The  subroutines  which  follow  were  written  in  ASCII  Fortran 
Level  9R1  (Sperry  Univac  Series  1100).  This  language 
allows  IF-THEN-ELSE  blocking  statements.  If  the  program 
is  to  be  used  with  a  version  of  Fortran  which  does  not 
permit  these  statements,  conventional  IF  statements  must 
be  substituted. 


Subroutine  ELPROP 


l 


i  .  .  . 


l  ’ 

i j  i  .... 

I  T 
i 

u 

i 5  c .... 

ib 


3 0  L 

7 1  r. 


34 


3  b 


THIS  ’  C'F’UU  ■  .  ^  DETERMINES  MATRIX  COEFF  1C  IF:! I  T 3  OF  L  LO  FOP 
I S0T"  ;:p  if  OP  OR  1H0TPCPIC  MATERIALS  UNDER  PLF^IL  STRESS  STPA  IN 
LOriiiiu  CONDITIONS  FROM  THE  hPFPPPIhTE  li.4;  "  !rL  U'NYThNTS. 

IfFLIC  IT  DOUBLE  PRECISION  ‘fi-H.O-O 
LOG  1  C  -  PPT 

r ot if  ; !  c dhTh  ■■  o , head ( 20 ) . nijmi ip ,  nui  fl  .  run  imot,  he r i .  1  iro .  j pp 

COMMON,  ELDATA  DMMMIh.MCT.  IEl  .  NEL 
com;  10: 1  Lip  I  IT  •  FRT 

[M : -F N S  1 0 H  I-  i.  1  *  . L ID  2  * , T >  3.3  )  .  D XL" 3 , 3 'j 
DMA  IJD  4HPESS ,  4HRA  IN 

READ  1 5. 10G0)  MATVPE- I .L.K .ANGLE 

READ  '  5  ■  100 1  D  •.  4 )  ..El.  PNU 1 2  .  E3  .  G .  E3,  PHU 1  3 .  PM  1.1 23 

DETERMINE  WHETHER  THE  PROBLEM  IS  FLhME  STRESS  OR  PLANE  STRAIN. 

IF  1  I .  HE  .  0  1  l  =  I 
IF  (  ]  .  E 0 . 0 >  I  - 2 

STOP  E  VALUE  G~  I  in  TYPE  »:FOP  USE  IN  SUBROUTINE  ELMTO 1  ■ 

D  (  1  ‘  =  I  lAT’VPE 

STORC  ORDERS  OF  GAUSIAN  QUADRATURE  TO  PE  USED  IN  STIFFNESS  MATP I 
AND  SI  PESO  EVALUATION  >L  AND  Y.  -  FOR  USE  IN  SUBROUTINE  EL.MTGP 
L  =  !  i  •  NO 1  3  ,  I  IAXm  1  1  .  L  1  1 
l  •  2  >  --  L 

L  -  MINO'  34  IP  MO-1  1  4  4» 

D  1  3  1  *  n: 

I  r  :  i-  r  T  .1  LIP  I TE  1  6 . 2000  *  IJD  y.  I ') Mm TYPE 

GO  TO  CURRENT  APPhY  PROCESSOR 
Lit  j  1  u  1  1 J ..  2 . 4  ■  .  I  in  I  YPE 

ISOTROPIC  MATERIAL 

14  5)  -  E 1 T ( 1 .  +  ll-P  f:RNU12M(l.  •**  PHU  12) /Cl.  -  ItRNU  12  1 
D  >: 6 j  *  FNU12*Di5Mf  1 .  +  >1  -  I)*RNU12) 
i'"3)  *  El/2,. ‘(I.  -}  RHU 12) 

IK?)  *  D ( 5 ) 

IF  1  ART)  WRITE  <S,  2101)  EURNU12 
GO  TO  20 


7  0  C 


40  C  .  .  , 

4 1  C  .  . 

42  2 

.  .  ORTHOTROPIC  MATERIAL 

I F  r  I  IhTYF'E  .  EQ .  2  .  AND .  PRT) 

WRITE  (6,2202 

43 

I F  ■:  M  A  TYPE .  EQ .  3 .  AND .  PRT) 

WRITE  (6,2302 

44 

IF  (PRT)  WRITE  (6.2203) 

E1,E2,RNU12,G 

45 

Dr 8)  -  G 

4b 

RNU21  =  RNU12*E2/E1 

43 

IF  (I.F0.2)  GO  TO  10 

48  C  .  .  , 

PLANE  STRESS 

49 

DUM  =  1.0  -  RNU 12*RNU2 1 

50 

DCS)  =  E1/DUI1 

51 

DCF)  =  E2'DUf1 

angle 
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I*  *3..i  =  PNui :  .  i*i  r  i 
GO  TO  20 

C -  PLANE  STRAIN 

10  RNU7 !  -  PN'.U  v‘F_7  E  I 

PMil  VO  -  «•  NU27  fr **’  *  fc- 2 

I'lJi!  FNi.i  !  2  t  h  N02  i  -  PM! !  1 ~  !  IU  3  1  -  Pt  l.?’*!' tPNNT?  - 

1  FN!  12  1  •  r !  !i  i’7.:  rPNi  1 1  7  -  F'NIJ  !  2  :  ?  NO  I  :  v'f  'MV  *1 

I’ »  3  1  -  F  :  »  •  1  .  -  P!!023>PilM73  ?  L*TN  1 

E"  b  ■  -  F.2  +  '  F'l  IN  1  2  *-  PNU72 < FNU 1 3 1  DNM 

I"  .  «  =  E2--  *  1  .  -  PNL!  1  3  *PN!  I  :<  1  ■  I'UI  I 

IF  t  F-T  .  UP  I  TF-  *0.2204'  L3 ..  F  NN  1 7  •  PMU23 

30  I F  f  PP  I  '  !.P  I  TF  1  0 . 200  1  •  I"  4  »  .  L  •  i_  .  f  i. 

IP  1  MAT  i  F'F  .  E'J  .  *■  1  ijij  T ! j  i< 

PE  T'J  PM 

C _  OP  "'HOT  POP  IC  l,J  T  rn  PE3PF.lT  TO  1-2  COORD  I  NATES 

C _  DETERMINE  TFANSFOPf ihTIOH  MATRIX 

3  C  ALL  TPhMSF  ■'  ANGLE,  T  ' 

C _  TRANSFORM  PROPERTIES  MATRIX  TO  GLOBAL  (>’,Y)  COORDINATES. 

CALL  BTPEP  1  T .  3 .  D  1  S  '  ,  D  £  1 ,  D < ?  1  -D<8'>  .  I' MY) 

i  -  _tr 

T '  0  31  1  =  1.3 
DO  31  J  = I  .  3 
D 1  k  )  =  d::Y‘.i,j"' 
i  =i+i 

31  CONTINUE 
RETURN 

r.  ,  .  .  . 

C . . . .  AN  I S OTP OP  I C  MATER  I AL 


4  IF  (PRT''  UPITE  »  6 . 2402 '«  El,  PNU12.E2.G.E3.PNU13 
D ( 5 )  -  El 
D i  6 .*  -  PNUI 2 

run  =  E2 
D(8)  -  G 
1 1 1 9 -  E  3 
D  '10)  =•■  PNUI  3 
GO  TO  20 

C _  FOPMAT  STATEMENTS 

' F0O  FORMAT  (4 15, F 10. 01 
1001  rOFMAT  18F1U.0) 

iCC’O  FOPMAT  I,'/.  51 '.8HP.ANE  ST.A4.23H  LINEAR  ELASTIC  ELEMENT,  10X.  '  MATYPE 
1=  •  .  12. 

2021  FOR!  'hT f  '20X.  7HDENS ITY.  101  :.E  18 .5,  •//,20X, 44HGAUSS  PTS  IN  XI,  ETA  PIPE 
lCTIutiS  RESPECTIVELY,,  .301-:. ’(I).  FOR  STIFFNESS  COMPUTATION*  »5X,  i  1 . 

2  * . ' .  II.  ’.30V. ’ (2) .  FOR  STRESS  COMPUTATION" .8X. 1 1 .  ‘  .  * »  1 1 ) 

2 i O i  FORMAT  i  10X. " ISOTROPIC  MATERIAL  * ,  , 1 5K, 1 8HY0UNG ‘ S  MODULUS  «  .  4X. 

1E18.5.  . 15: 1. ISHPOISSON’S  PATIO  =  . 10X.F8.5) 

2202  FORMAT  > 10X. " OPTHOTPOP IC  MATERIAL.  PRINCIPAL  DIRECTIONS  t  1.2'  COIN 
1C  IDE  1,1 1 TH  GLOBAL  lX.Y'  AXES  RESPECTIVELY*  > 

2203  FORMAT  UX. // . 15X.  THE  1  =  ,  E 18.5-/. 15X.7HE2  *  ,  E18.5./. 15X 

1 ,  THNU  1 2  =  ,  bX.  F8 . 5  •  '  •  1  OX,  TNG  12  *  ,E18.5> 

2204  FOPMAT  •  15X,  THE 3  =  .  E 3 8 . 5 ,  4  1 5X,  THNU 1 3  =  ,  6X,F0.5,  U  15X,7HNU23 

1  *  -6 Y , F  8 . 5 1 

2302  FORMAT  (  1 0 !•  i .  ’  0  P  TH  0  TP  OF*  I C  MATERIAL.  PRINCIPAL  DIRECTION  1  ORIENTED 
1  AT  A  CCU  ANGLE  OF  *,F9.4.‘  DEGREES  FROM  THE  GLOBAL  X  AXIS.’) 

2402  FOPMAT  ♦  10X.  ’ANISOTROPIC  MATERIAL'  .  15X.  'E'l,  1>  =  *  ,  E  18 . 5 .  U  15X 

1  *  '  E  *  1 . 2  ■'  =  UEIB.V  ,  15>:.  *E(1 .3'  =  *  .  E 1 8 . 5 .  U  5X .  '  E  <  2 . 2  •  =  ',E!8.5. 

2/.  15X.  'E'2.3i  -  '  .E18.5,/',  15X, -E13.3)  *  V.E18.5J 
END 
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84 
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♦  42 

♦  ro 
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90 


7  ■ :  o  i  -  ♦ 7 .  * 

\V*/:  ♦  li;  1 

7 "0"  i:  ♦10” 

jjr.i  ♦  ;  i  j‘  . 

7'  :>  7:  ’  lv  >  i  0  7 

7407  :  7  »I0'‘ 

in-  VARIABLES  *+* 


hML7  E 

f  i  7 

70 

[■■TREE 

7 

; 

1 0 

■*  14 

1 9 

+23 

♦  7  1 

*  cr 

♦  '  •'  '"-i 

♦  59 

♦  80 

+6 1 

Df  f 

Iuii 

♦  43 

50 

5  i 

.  C  “* 

59 

El 

1 0 

J  14 

73 

*;c 

37 

44 

E2 

>14 

44 

46 

51 

56 

E3 

♦  14. 

c«3 

5  b 

62 

83 

EL PROP 

1 

G 

+  14 

45 

0*7 

Og 

87 

HEAL 

■J7 

I 

+4  3 

*16 

*17 

26 

33 

IEL 

o 

;fp 

c 

.1 

+75 

7b 

K 

*13 

*24 

25 

63 

*73 

L 

*13 

*22 

■<7 

63 

MA 

8 

MhTYPE 

+  13 

19 

26 

•  29 

42 

MAX0 

22 

24 

MCI 

o 

MING 

~ 

24 

N 

8 

NEL 

8 

MEN 

7 

NECf 

7 

Hur  IEL 

7 

IiUMMhT 

7 

N  LI  MNP 

7 

U 

PRT 

( 

6 

9 

26 

37 

42 

PNU12 

*14 

33 

34 

35 

37 

PNU13 

+  14 

55 

57 

60 

61 

RNU2  1 

+46 

49 

57 

RNU23 

+  14 

56 

57 

59 

62 

RNU31 

+55 

57 

61 

RNIJ32 

*56 

57 

59 

60 

T 

10 

70 

72 

TRANSF 

70 

IJD 

10 

*1  I 

26 

+25 

+33 

*34 

*35 

*36 

*45 

*50 

63 

{  2 

+76 

*84 

*05 

*86 

*87 

60 

61 

46 

50 

cr  c 

j . j 

59 

83 

84 

60 

61 

83 

86 

og 


34  47  *74  75  76 


76  *'77 

43  64 


43 

44 

62 

63 

83 

44 

46 

49  . 

52 

57 

60 

83 

85 

62 

83 

89 
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Subroutine  TRANSF 


1 

SUBROUTINE 

TRANSF i BETA. T 

2 

IMPLICIT  DOUBLE  PRECIS  1 01 

7 

DIMENSION 

T '  3 , 3 ) 

DF 

<TA  FD  3. 

1 4 1 59265353079 

cr  r- 

THIS  SUBROUTINE  COHSTPljiI 

A  C  .  .  . 

FOR  AN  ANGLE  OF  BETA  DEG! 

7  c . . . 

orientated 

BETA  DEGREES 

3  L  . . . 

■V  COORD/ 

JhTES 

Q 

BETA  -  BETA+F'I  1B0. 

10 

T  i 

i .  i  :■  =  d 

I  OS  (BETA) 

11 

T  i 

1 .  ?  i  =  D 

SIN  (BETA) 

1 2 

T> 

1-3)  -  D 

SIN  (BETA)  *I«COS 

13 

T  i 

2,1)  =  7 

«:  1,2) 

14 

T  i 

2  «  2  1  -  T 

'  l ,  l ) 

15 

T< 

2 . 0  kT  i.  1 , 3 ) 

1  b 

T 

2  •  3  “  — 

1 .0*T k 1 , 3) 

17 

T 

3,1)  = 

2 . 0-f-T 1  2,3) 

18 

T> 

3,3)  =  T 

a,n  -  T < i , 2 ) 

1 9 

RE 

:turn 

20 

END 

BETA 

i 

*9  10  11 

DC  OS 

10 

12 

DS  IN 

!  1 

12 

PI 

9 

T 

1 

3  *10  *11 

TRANSF 

1 

N  ».h-H . 0-Z'1 
3238D000  ' 

TS  THE  PLANE  TRANSFORMATION  MATRIX  CT3 
PE EC .  I.E.  THE  LOCAL  1-3  COORDINATES  ARE 
COUNTER  “CLljr L LI  1 SE  OF  THE  GLOBAL 


VARIABLES 

12 


i 

L 


Subroutine  BTREB 


SUBROUTINE  BTREB  (B,  1BCOLS.E  M -E  12.  E22.  ET'i.F'POD') 

THIS  SUBROUTINE  COMPUTES  THE  MATRIX  MULTIPLICATION  BT+'EyB  =  PE 
E  IS  A  3  BY  3  SYMMETRIC  MATRIX  WITH  HCH  ZEPU  COMPONENTS  E(1.U 
F  1  1 . 2  ' .  E  ( 2  •  l1 .  E  1 2  >  2  V.  E  <  3  •  3  < .  B  IS  A  3  Fi Y  IBcuLG  MATRIX. 

PFCiTi  IS  A  SYMMETRIC  MATRIX.  ONLY  THE  UPPER  MALI'  IS  COMPUTED. 


.  .  TH I  'I 

SUBROUTINE 

: 

..El  7 

h  7  BV  3  b 

j  r .  . 

.  .  F  '  1  . 

2>  .E(2. 1 ’ .E 

L* 

.  .  F-’F ' Hi 

IS  A  SYMME 

f. 

I  r  ifl 

IC IT  DOUBLE 

1  ' 

D  J  ME' 

IS  ION  B'3. 

-• 

FO  1 

1  =  1.  I POOLS 

r« 

B  1  I 

=  Elf  1  .  I .' 

\  >  ‘ 

tv:  i 

--  E: '  2 . 1' 

*  1 

Burn 

*  E11+B1I 

,  “ 

"DU  r  12 

*  E 12*911 

17. 

Dl.il  17 

*  E33+EC3, 

14 

DO  1 

J=I . 1BCOLS 

15  1 

F'POD 

<  1 ,  J  <  =  Em  1 

RE  I  URN 


1 

o 

14  +15 

El 

1 

7  9 

B 1 1 

¥3 

1 1  12 

B2I 

MO 

I  1  12 

BTREB 

1 

BUM  1 

M  1 

15 

DU  M2 

+  12 

15 

BUMS 

+  13 

15 

Ell 

J 

1 1 

E 12 

1 

li  12 

E22 

1 

12 

E33 

1 

13 

I 

+8 

9  10 

IBCOLS 

1 

7  S 

*•**  STATEMENT  NUMBERS  *** 

***  VARIABLES  *** 

?  10  13  15 


SubfOb  ^\e  SHAP3 

••  'j.E.p  ;  r.H-’P'  !  "  P:- .  ^ .  i: :  ’ 

,■  1  , .  :  t*  i_l  t=  U  t :  T I  NT-  i.ll  ‘  •  i'l  r  ■ !-  l  F  i  -  n’  T  I  f.ii'-j  AND  N77E3SAPY 

• .  ■ ■  ;v.iTi\r' .  f..r  m:  : , / 1 n- r'  :c  \ub  r-:x  '--m  peocped 

■  *  ■  •  {_  T  ]  ]  1  1  1 i  1 1  i  r ;  \  ?•  j-  ;•  r  :  |— i  —  H'  k  i  t  —  , 1 

;■  .....  v  r  r  ‘  .'b.- {_  f  >r  ~  c  n t*  f  «f p r  c  r .  r  ;j t ?  -■  -p^f :  y 


:  “  ■  1  .  -  E  Th  *  ETA  '  I' . 

i0  qa  -  :r.  it. 

n  QUETH'  -  1.  3. 

:r  c  :cn  -  cc;s:: 

17  C  FTA  r  7  =  C  *  J  » T  7; 

c  i;  ;  -  -  onethd 

17  =  XI  +  ONETHD 

17  ClE'i'H  -  ETh  -  ONETHD 

1  r  C.7l:  TA  a  EM  +  CNETHD 

iS  DLN  =  C 1 1 1  SO  -  COrXl*CiXI 

1 9  Ij 7;: .  I  t  C : •!  1 5 0  -  C  0  + X I +-C 2 7! I 

20  Tit*:-  «  CETASO  -  COWETA iCIETA 

21  [7Eth  -  C  ETH  SC  -  CG  +ETA4C2ETA 

22  C  ININ  =  1.  -  XI 

23  C  « r : :  i  -  i .  +  xi 

2-i  CU1  ETh  «  1.  -  ETh 

?c  Cl P ETh  ■-  1,  +  ETA 


■m  : _  ZERO  SHF  >  J,  I  lit 

22  DO  I  HQ  I  *  5,12 

72  DO  i  OO  J  35  1.-3 

50  ICO  SHPiJ-I)  *  0.0 

51  C _  IF  D IDS  IDE  NODE  II  IS  NOT  SPECIFIED,  I.E.  IXCII)  =  0, 

32  C _  THEN  SHP'J.Il)  ,J  =  l-3,tJJLL  STILL  EQUAL  0  AFTER  THE 

SS  C -  FOLLOWING  STATEMENTS 

3-  C _  THUS ,  THE  CORNER  NODES  ON  THIS  EDGE  WILL  RETAIN  ONLY 

3?  C -  QUADRATIC  OR  LINEAR  TERMS 


L uHFUTE  SHAPE  FUNCTIONS  AND  DERIVATIVES  -  SEE  FIG. 2  FOR  NUMBERING 

IF  ( 1 1 1 1  5 > . EQ . 0 . AND . IX( 9) . EQ . 0)  GO  TO  101 
IF  ( I K *  5 > .HE. G. AND. IX(9) .NE.dV  THEN 
SHF i;  I.5-  =  -C  inETMDIXl 
SHP (2.5)  *  CiXI*CXISQ 
SHP  ( 3 , 5 )  -  -Cl ME TA*SHP  (2,5) 

SHP (1,91  *  C 1 META *B 2X1 
SHP  (2,9 1’  *  -C2XI*CXISQ 
SHP  v  3,9;  23  -Cl  META  49  HP  ( 2,91 
ELSE 

SHP  f  i .  5  >  =  -X I  <•  C  i  META 
SHP j12.5>  *  -32 
SHP  (3.5 )  *  S2;<  7  1!  iETA 
END  IF 

IF  t  I X ( 61 . EQ . 0 . AND . I X ( 1 0 ) . EQ . 0 1  GO  TO  102 
)r-  ’1  I X ( 6 1  . HE  . 0 .  AND  .  IX'  10)  .NE.C)  THEN 
31  IP '  1 . 6 >  =  -C  IE  f A  *C  ETASQ 
SH'P  k  2,  bi  =  -C1PXMD1ETA 
SHP ‘3,6;  -*  C 1 PX I  ^  SHP  (1,6) 

S  HF*  '1.1 0 1  =  C  2  E:  T  A  v  CETh  S  Q 
SH° K 2 ,  1C'«  *  C1PXUD2ETA 
SHP'.'SmCt  -  CliXMSHF’U,  10> 


9 


b4 

65 

66  102 


HP  ‘ 

1 . 

6.' 

=  T2 

HP 

v . 

6 ' 

-  -ETA*C1P 

XI 

■  HPC 

3  • 

6  ‘ 

=  T2+C 1PXI 

IF 

7  i 

EO 

.0. 

AND. IX' 1 1 ' 

.EO, 

,  n  i  go 

IX< 

7' 

.HE 

.0. AND. I X ( 

IP. 

.  NE . 0 • 

Cl. 

X 

1, 

7  i 

=  ClPETAfD 

2X1 

HP  C 

-i 

7 ;» 

-  C2XI  ♦C<‘I 

SO 

HP ' 

. 

7  ,i 

-  CIPETA+S 

hp  ' : 

HP'. 

1, 

1  P 

=  -C 1FETA 

rVl: 

:i 

HP', 

2 . 

IP 

=  -cixpr 

XiSC 

j 

HP  ( 

3, 

1  P 

-  C IPETh  * 

SHP- 

2 ,  1  1  '« 

HP  i. 

1 . 

7) 

-  -xpcipeta 

HP  < 

2 . 

7 j 

=  S2 

HP  ( 

T  _ 

“•  *) 

=  S2h:ifet 

A 

IF 

8  1 

EO 

.0  . 

AND. IX i 12' 

.EO, 

,  O  -  GO 

IF  'IX (8 
SHP  el 
SHP  (2 


yu 

91 

92  C . . . 

93  C . . . 

94  104 

95 
9b 

97 

98 

99  C  .  .  . 
100 

101 

102  111 

103 

104 

105  121 

106 

107  C. . . 

108 
109 

11G  112 
1 1 1 
1 12 

113  122 

1 14 

115 

116  109 
1  17 

1  18 


.  HE  .  0  .  AMD .  IX«  12  '  .  HE  .  0  '  THEM 
8  >  =  -C  2E  Tfi  * C  ETAS  0 
8)  *  c  imxi  +D2ETA 


SHP (3. S' 

=  -C 1 MX I + SHP 'IX 

SHP '  P  12 

C  1ETA  mIETASO 

SHP- 2,  12 

-  -C1MXPD1ETA 

SHP (3.12: 
E 

=  -Cl MX I  *SHP(.  1 

SHP* 1-3) 

=  -T2 

SHP i 2.3) 

=  -ETA+C  1MX1 

SHP (3. 8) 

=  T2*C1MXI 

END  IF 

CORRECT  CORNER  NODES  FOP  PRESENCE  OF  MIDSIDE  NODES 
K  =  8 
KK  =  12- 

DO  109  I  =  04 
L  -  I  +  4 
LL  *  I  +  8 

ADJUST  CORNER  SHAPE  FUNCTION  FOR  NODES  ON  THE  CU  SIDE 
IF f IX(KK) .HE. 0. AND. IXOO  .NE.0)  THEN 

DO  111  J  =1-3 

SHP(J.I)  =  SHPfJ.P  -  0NETHD*l2 . 0*SHF t  J  #  KK)  +  SHP(J.K)) 
ELSE 

DO  121  J - 1 . 3 

SHP ( J >  I  '  -  SHPi  J, I '  -  0. S^SHPr  j.Kt 
END  IF 

ADJUST  CORNER  SHAPE  FUNCTION  FOP  NODES  ON  THE  COJ  SIDE 
IFC IX(LL) .NE.0.AMP. IX»L  » .HE. 0»  THEN 

DO  112  J  =  1  *  3 

SHPCJ.H  *  SHP  i  J .  I  1  -  ONE  THD ♦ ' 2 . 0  +S.HF  U .  L  1  +  SHP(J.LL)) 
ELSE 

DO  122  J-l .3 

SHP(J.I)  =  SHPiJ.  P  -  O  .  v  ♦SHF  «  J  .  L %t 
END  IF 
K  =  L 
KK  =  LL 
RETURN 
END 


f  ♦  t 
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100 

j‘b 

>30 

101 

40 

*57 

102 

Cj-? 

+  tb 

1 03 

t'D 

+  73 

104 

79 

f  4 

109 

9  b 

+4  lb 

1  t  1 

101 

■+  102 

1  12 

1 09 

+  1  1  0 

12  1 

104 

+  105 

122 

1  12 

+  113 

*  *:•-+ 

VARIABLES 

•f:  f  +: 

CO 

+  10 

12 

13 

13 

19 

20 

?1 

ClETA 

1 1 6 

2+3 

crer 

34 

C  1MET 

+  24 

42 

44 

45 

47 

49 

51 

cm::: 

+  22 

r!  3 

O  T 

C'C 

86 

89 

90 

C 1PETA 

*  2 1 

68 

70 

71 

77 

7*  S 

(  l 

C  ip:  A 

+  2" 

5  b 

c  — 

59 

60 

63 

64 

c  i : :  i 

+  1 4 

IS 

43 

72 

C2ETh 

+  1 7 

21 

n*  o 

31 

C2: :  I 

+  15 

19 

46 

69 

C ETA SO 

+  1 7 

20 

21 

c:  Cj 

rr  r 

81 

34 

C  ■  •  I  so 

+  12 

18 

19 

43 

46 

69 

7  "■ 

PIETh 

+  20 

5  b 

O  cr 

pi:  :i 

+  13 

42 

7  1 

P2ETA 

+  21 

59 

82 

P2: :: 

+  19 

45 

68 

ETA 

t 

9 

16 

17 

20 

21 

24 

*?cr 

63 

89 

i 

+  23 

3G 

*96 

97 

93 

102 

105 

130 

1 13 

i: : 

1 

5 

40 

41 

c  *7 

54 

66 

67 

79 

80 

j 

+  29 

30 

*101 

102 

*404 

105 

r  [  09 

i  10 

*4  12 

113 

K 

+  94 

100 

102 

105 

*115 

1  1 

*95 

100 

102 

+4  16 

L 

+  97 

IDS 

1  10 

1  13 

115 

LL 

+98 

108 

1  10 

1  1 6 

GMETHP 

+  1  1 

14 

15 

1 6 

17 

1+02 

ne 

+•  3 

12 

50 

51 

7  b 

77 

SHAP3 

1 

SHP 

1 

tr 

+  30 

*42 

*43 

*44 

*45 

*46 

*47 

*49 

*55 

*56 

+  57 

*58 

*59 

*60 

*62 

*63 

^64 

*68 

*71 

*72 

*73 

+  75 

’  *  7  6 

+  77 

*61 

*82 

*83 

>-o4 

+  33 

*89 

*90 

*402 

*105 

*1  10 

*4  13 

T2 

*9' 

13 

62 

64 

88 

90 

XI 

1 

8 

14 

15 

18 

19 

22 

23 

49 

-?ct 

100 


*50 
*6  9 

-C;  c: 


108 


*5 1 
*70 
*86 
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Subroutine  ELMT01 

1  SUBROUTINE  ELMTOl'D.ULXL,  IX, TL -S, P, HDF . MDM. NST.  ISLJ' 

2  C 

2  C _  FLAME  STRESS- 3 TP A  IN  ELASTIC  LINEAR  ELEMENT  PONTINE. 

4  r 

5  IMPLICIT  DOUBLE  PRECISION  iA-H.O-Z' 

6  COMMON  CDhTh  0. HEAD'  20  ■  .IIUMNP .  NUMEL  .  NIJMMAT.  NEN .  NED .  IFF' 

Z  COMMON  'EIDATA  DM. H . MA . Ml T. IEL . MEL 

DIMENSION  D  -  1  '  JJL'NDF.  I  '  X’L'NDM.  !  •  .  IX'  1  »  .  TL  ■  !  :  -S'hST.  1  •  .PCI  » 
0  l  .  SHP  i  3  .  1 2  ' ,  X  I  i  9  i  .  E Th  1  9  1  .  M  i  »  9  •  .  S  1 0  ■’  0  »  .  E  F  :  1  “  ► 

:  0  C  .  .  .  . 

11  C _  Gu  TO  CORRECT  hFPhY  PFOCESSOP . 

12  SO- TO  ‘  1 .2.3.  4.5.4' .  ISU 

i  t  r 

14  C . . . .  INPUT  MATERIAL  PROPERTIES. 

15  C _ 

It  I  CALL  ELPPOP'D  ■ 

IT  LINT  =  0 

1  P  E  !  U  F1  M 

10  2  Rt TURN 


IUE  ELEMENT  STIFFNESS  MATPIX  SCI,,!' 


31 


40 

41 


L  =  D 1  2  ‘ 

IF  r  LML  .  NE  .  L  !NT  ■  C^LL  PSPUSS  'L  .  LINT,  XI  .  ETA,  LLP 
MhT'.FE  *  D  ‘  I  * 

POP  EhCH  INTEGRATION  POINT  COMPUTE  STIFFNESS  CONTRIBUTION 
r«0  320  L  =  I.  LI  NT 

CAtaL  SHAPE  '  "I '  L • .ETA' L ■ . XL . SHP, DETJAC , NDM, NEL . IX, .FALSE. ) 
E;'v  =  I'ETPAC  NJT ( L )  ■ 

[•  1  1  -  D  ■  5  :  +T'V 
D 12  =  D'6'fDV 
IF  M  -A TYPE  . LE . 2)  THEN 
1*22  =  D  >  Z  ’<  h  D V 
D33  -  DcSVtDV 
ELSE 

D  1 3  =  D (7)*DV 
D22  =  D>8)*DV 
D23  =  D»9)*DV 
D33  =  DU0)*DV 
END  IF 


FOP  EACH  NODE  J  COMPUTE  DB  =  B*B 


44 
.  fr“ 
4b 

M  ; 

4;S 


50 

C.  1 


53 

54 


SB 

59 

60  C. . . . 

61  C . .  .  . 


DO  320  J  =  l . NEL 
IF  < MATYPE.LE.2i  THEN 
DB 1 1  *  D 1 1 *SHP 1 1 . J » 
DB  12  -  D 12  f  SHP  <"2  ,  J  ’■ 
DB2  I  -  D 12+SHP  <  1 ,  J) 
DB22  =  D221SHPC2. J) 
DBS  1  -  B 3 3* SHP (2 . J) 
DB32  =  D33  +  SHP  U ,  -O 
ELSE 

DB  1 1  =  D 1 1  -fSHP  ( 1 ,  J  ' 
DB 12  *  D 12+SHP  C 2 -  JO 
DB2 1  =  D 1  2*SHP  r  1 .  J ) 
DB22  =  D22*SHPC2,J> 
DB3 1  =  D33  *  SHP ( 2 , J ) 
DB32  =  D33*SHP  U ,  J) 
END  IF 


+  D  13*5HP(2,  J) 
4-  D13*SHPC1,J) 
+  D23*SHP<2, J) 
+  D23*SHP< 1, J) 
+  B13*SHPC1,J) 
+  D23*SHPC2, J) 


FOP  EACH  NODE  I  COMPUTE  S  *  BT*D8 


12 


-  Hf  " 


r 


b3  32m 
63  C . . . 

ro  c . . . 
r  i  c . . . 


L  .  .  . 


DO  32D  1  -  1  .  J 

S»  I  +  I-l  .  J+-J -  1 .*  =•?.«■  ]-!*:-]  .J+J-l  *+SHP'  1.1'  4  DB1 1+SHp!' 2.  I  ■'  *DB3 1 
3  1  I  f  I  -  l  .  .1  +  I  I  -  i  I  ■+*  I  "  1  .  J 4j  1  +'jnr  1  1  .  i  •  i!P  .i+bHP  1  J  •  I  1  *  PB34 
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i  15 

MCT  =  50 

1  lb 

430 
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CONTINUE 
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A  finite  element  program  is  presented  for  stress  analysis  of 
wood  members  of  arbitrary  shape  which  are  subjected  to 
plane  strain/stress-loading  conditions.  This  program  extends 
one  which  is  available  in  the  literature.  User  instructions  and 
a  listing  of  the  developed  subroutines  are  presented. 

Keywords:  Finite  element  analysis,  computer  program, 
isoparametric  elements,  stress  analysis,  orthotropic 
materials,  anisotropic  materials,  plane  loading,  design,  cubic 
finite  element,  quadratic  finite  element. 


